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A general relation is derived for the action difference between two fixed points and a phase space 
area bounded by the irreducible component of a heteroclinic tangle. The determination of this 
area can require accurate calculation of heteroclinic orbits, which are important in a wide range of 
dynamical system problems. For very strongly chaotic systems initial deviations from a true orbit 
are magnified by a large exponential rate making direct computational methods fail quickly. Here, 
a method is developed that avoids direct calculation of the orbit by making use of the well-known 
stability property of the invariant unstable and stable manifolds. Under an area-preserving map, 
this property assures that any initial deviation from the stable (unstable) manifold collapses onto 
themselves under inverse (forward) iterations of the map. Using a set of judiciously chosen auxiliary 
points on the manifolds, long orbit segments can be calculated using the stable and unstable manifold 
intersections of the heteroclinic (homoclinic) tangle. Detailed calculations using the example of the 
kicked rotor are provided along with verification of the relation between action differences. The 
loop structure of the heteroclinic tangle is necessarily quite different from that of the turnstile for a 
homoclinic tangle, its analogous partner. 

PACS numbers: 45.10.-b, 05.45.Ac, 05.45.Pq 


I. INTRODUCTION 

There are significant classes of problems in physics and 
chemistry for which it is very important to have accu¬ 
rate information about periodic, heteroclinic, and/or ho¬ 
moclinic trajectories in chaotic dynamical systems. For 
example, they arise in semiclassical trace formulae [1] 
and wave packet propagation [2, 3]. Trace formulae 
have been invoked to understand shell structure in nu¬ 
clei, metallic clusters, and quantum dots [4], and the 
Bohigas-Gianonni-Schmit conjecture relating chaotic sys¬ 
tems quantum spectra to random matrix theory [5] . The 
time evolution of wave packets have been used for un¬ 
derstanding driven cold atoms [6], electrons in strong 
fields [7], fidelity studies [8, 9], and a broad range of 
spectroscopic and pump-probe experiments [10, 11]. If a 
particular system of interest possesses a strongly chaotic 
dynamics, then a semiclassical approximation leads to 
sums over finite-time segments of such trajectories [3]. 
Furthermore, the quantum mechanical phases are largely 
controlled by Hamilton’s Principle Function (or classical 
action for short) for these segments divided by Planck’s 
constant. As a result, deep in a semiclassical regime, 
small changes in classical actions or small action differ¬ 
ences between the various trajectory segments may result 
in significant changes in interferences quantum mechan¬ 
ically. The actions must be known accurately to predict 
these interferences correctly. 

As time increases in the semiclassical sums, the num¬ 
ber of contributing terms increases exponentially rapidly, 
but there cannot be an exponentially increasing amount 
of information present in the quantum propagation. This 
situation has to be reflected in the classical dynam¬ 
ics through the existence of correlations in classical ac¬ 
tions [12]. A critical element that bears on the correla¬ 
tions is the relationship between the limiting differences 


of certain heteroclinic trajectory pairs and areas enclosed 
in phase space [13]. These areas show up in multiple ac¬ 
tion differences and their boundaries are determined by 
the heteroclinic or homoclinic tangles themselves. Our 
purpose is to relate periodic orbit or fixed point action 
differences to phase space areas defined by heteroclinic 
tangles, and develop a scheme in which heteroclinic or¬ 
bits of strongly chaotic systems are calculated accurately 
enough and over a long enough time interval that the lim¬ 
iting classical action differences give precise evaluations 
of such critical areas. 

It turns out that for a chaotic system with a large 
enough Lyapunov exponent, using Hamilton’s equations 
with any phase point on some particular orbit (or if it is a 
mapping that is of interest, the map) will fail to faithfully 
follow a heteroclinic orbit segment, except for very short 
segments. In order to do better than this, an alterna¬ 
tive computational scheme is desirable. Previous papers 
on this topic were mainly focused on continuous time 
systems, where the infinite time interval associated with 
the homoclinic or heteroclinic orbits are truncated into 
finite time domains, with appropriate boundary value 
conditions so that the standard boundary value prob¬ 
lem solvers apply [14-17]. Modified approaches using 
the arclength of the orbit as the system parameter in¬ 
stead of time were given in [18] and [19], which avoid the 
truncation process but still require solving the boundary 
value problem using numerical discretization and colloca¬ 
tion techniques. Methods using spectral expansions that 
also avoid the finite truncation were given in [20]. Re¬ 
cent papers facilitating other numerical techniques can 
be found in [21], where Hermite-Fourier expansions were 
used to approximate the homoclinic solutions; and in [22], 
where a variational approach similar to the one proposed 
in [23] was employed. The method developed here is 
vastly simpler than those mentioned above. It focuses 
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on Hamiltonian systems with invariant manifolds, and 
relies on the fact that a heteroclinic (homoclinic) orbit 
lies at all times at the intersection of a stable and un¬ 
stable manifold. These manifolds can be computed in 
quite stable ways, and therefore their intersections can 
also be calculated in just as stable ways. It is concep¬ 
tually straightforward, also extremely fast in terms of 
calculation time. No special numerical techniques such 
as collocations in the boundary value problems or alge¬ 
braic equation solvers are needed. Therefore, the degree 
of precision is only affected by the density of data points 
interpolating the manifolds. If the set of data points are 
dense enough, which can be easily arranged by insert¬ 
ing more initial points for generating the manifolds, the 
precision is only limited by the interpolation technique 
used to intersect the stable and unstable manifolds, and 
machine precision, whichever is greater. 

By switching initial intersections (denoted by Ro 
ahead), our method can be used to find any heteroclinic 
(homoclinic) orbit with arbitrarily long excursion length, 
not being restricted to some particular orbit. Thus, an 
enumeration of the infinite set of orbits is made possible 
by switching the initial intersections, i? 0 , and repeated 
use of this method. In general, there is a minimal or 
irreducible structure from which to select an Rq in or¬ 
der to obtain any orbit of one’s choosing and to avoid 
duplications. For the homoclinic case, this is the well 
known turnstile (assuming no special symmetries), and 
for the heteroclinic case a necessarily different loop struc¬ 
ture. The relation between phase space areas and peri¬ 
odic orbit action differences reinforces this. Without loss 
of generality, from here on, we concentrate on dynamical 
maps as a continuous dynamical system can be reduced 
to a Poincare map. Thus, a heteroclinic orbit segment is 
constructed as a sequence of manifold intersections for¬ 
ward or inverse in time (iteration number) from some 
given intersection phase point as opposed to the forward 
or inverse mapping of that same phase point. 

This paper is organized as follows, the next section cov¬ 
ers necessary background and the introduction of some 
useful notation. All of the main ideas are illustrated us¬ 
ing the kicked rotor with a strong kicking strength, a 
well-known, simple paradigm of strongly chaotic dynam¬ 
ics. The following section describes the use of the well 
known structural stability inherent in chaotic dynamical 
systems [24] in order to construct stable and unstable 
manifolds. This is followed by a technique to locate in¬ 
tersecting points. The technique is used to compare sta¬ 
bility exponents of certain fixed points of the map cal¬ 
culated using the stability matrix and the constructed 
orbit segments, and to compare certain phase space ar¬ 
eas with limiting classical action differences. The area of 
the heteroclinic loop structure is shown to be the action 
difference of the fixed points. Finally, we summarize and 
note some ideas of interest for future work. 


II. BACKGROUND 

Let X = (q,p) represent points in phase space and 
classical transport be described by a time-dependent cor¬ 
relation function between an initially localized Gaussian 
density of phase space points p a centered at point X a 
and a final destination Gaussian density pp centered at 
point Xp. In the limit of highly localized densities, it 
can be expressed as a sum over all heteroclinic trajectory 
segments [3] 

Tp a (t)—>Y,{p^Pc) (!) 
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where T is a linearized dynamical mapping of points X , 
t is the time, i.e. number of iterations of the map, and 
7 denotes those segments that take time t to leave the 
neighborhood of X a and arrive in the neighborhood of 
Xp (they serve as the orbit segments about which the 
linearizations are done). For simplicity, let X a and Xp 
be fixed points of the map T. The segments belong to 
heteroclinic orbits that converge to X a for t —» — oo and 
converge to Xp for t —> +oo. The quantum mechanical 
analog in the semiclassical limit of the correlation func¬ 
tion is [3] 

C/9a(t)«X> \U^t)\a) (2) 
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where |a) is the ket vector corresponding to a quan¬ 
tum wave packet centered at X a and U 1 {t) is an ap¬ 
propriately linearized unitary time translation operator. 
The summation is over the same heteroclinic orbit seg¬ 
ments. These two equations hold equally well for open or 
bounded systems. The quantum expression is our main 
motivation for investigating the heteroclinic orbits and 
their properties. 

For a dynamical map a complete orbit {i?o} 
can be represented by a bi-infinite sequence of the 
form: {R_ oa ,..., R_ 2 , R-i, Rq, Ri, R- 2 , ■ ■ ■ , R^} where 
Rj maps to the point on the orbit Rj+t after t iterations 
of the map. The times have a translational arbitrariness 
to them, but here they are set such that 0 is the time for 
a presumed known initial condition that defines a partic¬ 
ular orbit. Given that in calculations the mapping itself 
cannot be applied exactly and repeated mappings lead 
to exponentially growing errors, a true orbit cannot be 
found this way. Ahead, we will choose a case for illus¬ 
tration whose dynamics are so unstable (large Lyapunov 
exponent) that it is not possible to follow a heteroclinic 
orbit segment more than ±5 iterations via the mapping 
(using double precision). Because our interest is in the 
set of heteroclinic orbits, it is not possible to rely on the 
Shadowing Lemma. While there is a true orbit of the 
system ’’shadowing” an orbit found with the mapping, 
the likelihood that it is the actual heteroclinic orbit of 
interest is vanishingly small. 
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A. The kicked rotor 


The kicked rotor on a torus has been a simple, yet 
extraordinary paradigm for chaotic systems for roughly 
50 years and a great deal is known about it [25]. It is a 
mechanical-type particle constrained to move on a ring 
that is kicked instantaneously every multiple of a unit 
time, t = n. The Hamiltonian takes the form 

2 oo 

H{q,p) = y - ^2 cos2 7r<? S(t — n) (3) 

TI— — 00 

The mapping equations are: 

K 

Vn+\ = Pn - — Sin 2nq n (mod 1) _ 

q n + 1 = q n + Pn +1 (mod 1) 

As the kicking strength parameter K increases away from 
zero, the system becomes more and more chaotic. For K 
values exceeding roughly 27 t, the system is very nearly 
completely and strongly chaotic. The Lyapunov expo¬ 
nent A is known analytically to be [26] 

A ~ ln 4 + 77y^ + o(- 1 - A (5) 

2 A 2 - 4 \(K 2 -4) 3 J 

In all the calculations ahead, the kicking strength is taken 
as K = 8.25, which is a strongly chaotic case. In addition, 
there are two very convenient fixed points of the mapping, 
(0, 0) and (0.5, 0), to be used for X a and Xp respectively. 
Therefore, the heteroclinic orbits discussed ahead lie at 
intersections of the unstable manifold of the phase point 
(0, 0) with the stable manifold of the phase point (0.5,0). 

It turns out to be convenient to use the “unfolded 
torus”, meaning that by not invoking the modulus 1 op¬ 
erations in the mapping equations, there is a “flat” phase 
space that extends to infinity. Each unit square is a rep¬ 
etition of the fundamental torus which is the [0,1) x [0,1) 
square in the phase space. Any two points that are sepa¬ 
rated by integer numbers on either the q or p coordinates 
are the same point; i.e. (q,p) and (q + n q ,p + n p ) are the 
same point if both n q and n p are integers. The integers 
n q and n p can be thought of as winding numbers (includ¬ 
ing negative integers), i.e. how many times a particle has 
wrapped around the cycles of the torus, which can have 
phase consequences in quantum mechanics. 

As shown in Fig. 1 , we define the part of the unstable 
manifold (solid curves) of X a which is initially pointing 
to the upper-right direction to be the upper branch of the 
manifold, denoting it by C/ + (0, 0). Denote the piece of 
the unstable manifold of X a which is initially directing to 
the lower-left to be the lower branch U~( 0, 0). The total 
unstable manifold including the two branches is denoted 
by: (7(0,0) = U + (0, 0) (J U~ (0,0). Similarly, denote the 
part of the stable manifold (dashed curves) of Xp which 
is initially directing to the upper-left to be the upper 
branch 5 + (0.5,0). The lower branch A - (0.5,0) and the 



FIG. 1. The initial segments of the unstable manifolds (shown 
in solid curves) of (0,0) and (1,0), and the stable manifolds 
(shown in dashed curves) of (0.5,0) and (1.5,0) on the unfolded 
torus. Ro is a heteroclinic intersection point between U + (0, 0 ) 
and S + (0.5,0). For all the figures in this paper, the unstable 
manifold is plotted as solid curves and the stable manifolds 
in dashed curves. 

total stable manifold 5(0.5,0) = 5 + (0.5, 0 ) (J 5“ (0.5, 0 ) 
are denoted in the same way. 

The motivation for distinguishing the upper and lower 
branches stems from the reflective property of X a . Each 
iteration of Eq. (4) will map t/ + (0,0) into U~ (0,0), and 
vice versa. Thus points on 17(0,0) “jump” between the 
two branches with iterations. Define the twice iterated 
map to be the compound map of two successive map¬ 
pings under Eq. (4), then t/ + (0,0) and U~ (0,0) are in¬ 
variant. On the other hand, Xp is non-reflective, points 
stay on the same branch of the stable manifold with it¬ 
erations. As an example, the heteroclinic orbit segment 
starting from R a (R 0 G U + (0, 0) fj 5 + (0.5, 0)) in Fig. 1 
will remain on U + (0, 0) (~) 5 + (0.5,0) for even iterations 
R± 2 , R± 4 , •••, etc. For odd iterations it will jump to 
U~(0, 0) p| 5 + (0.5,0). In the following sections the or¬ 
bit segment {Rj : Rj+t} is considered and a method is 
developed to obtain large numbers of iterations in both 
forward and inverse directions of the mapping, which is 
inaccessible via straight calculations of Eq. (4) because 
of the exponential growth rate of initial error. 


III. THEORY AND CALCULATIONS 

A. Stability in the neighborhood of invariant 
manifolds 

It is well-known that regions near stable and unsta¬ 
ble manifolds inherit strong stability in the sense that 
any point in the neighborhood of the stable (unstable) 
manifold with a small deviation from the manifold will 
exponentially collapse toward the manifold under inverse 
(forward) iteration [27]. It guarantees the legitimacy of 
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many methods to generate the global stable and unstable 
manifolds by propagating points from the corresponding 
subspace near the fixed points. The invariant manifolds 
can thus be calculated to high numerical accuracy and 
in quite stable ways, as shown by earlier researchers fa¬ 
cilitating different approaches [28-42]. A comprehensive 
survey of methods can be found in [43] . Consider the ex¬ 
ample schematically illustrated in Fig. 2, where a portion 
of a general homoclinic tangle commonly seen in Hamil¬ 
tonian systems, formed by the unstable manifold of a 
non-reflective fixed point Xp and its stable manifold is 
illustrated. Note that Fig. 2 is not specific to the kicked 
rotor. Heteroclinic tangles for the kicked rotor are sig¬ 
nificantly different from Fig. 2. However, the structural 
stability of the manifolds is always maintained. Suppose 
the point Xq is a true point on the manifold and X is 
found as a result of some numerical calculation instead of 
Xq . The deviation of X from the stable manifold shrinks 
to zero under inverse iteration, making X collapse ex¬ 
ponentially toward the stable manifold. However, the 
tangential deviation between X and Xq along the stable 
manifold will be magnified exponentially, and the two 
points will be spaced further apart along the manifold. 



FIG. 2. Manifold stability: deviation from the stable man¬ 
ifold collapses toward the manifold under inverse iteration, 
respectively. Assume for illustration that the phase point Xp 
is shown wrapped once around the torus, and is the same 
point in both locations. 


To show this, imagine performing a normal form trans¬ 
formation of the heteroclinic tangle on the left in Fig. 2. 
The picture in the normal form coordinates is shown in 
Fig. 3. The P and Q axes are the images of the unsta¬ 
ble manifold and the stable manifold of Xp respectively 
under the normal form transformation. The convergence 
zone of the normal form series extend along the mani¬ 
folds to infinity [44]. Assuming that A' lies inside this 
convergence zone, it evolves under the equation [44, 45] , 



FIG. 3. The image of the stable and unstable manifolds of 
point Xp in the normal form coordinates. The invariant hy¬ 
perbola under the mapping of Eq. (6) that passes through the 
point X is also plotted in dashed curves. 


U(QP) 

Q' = U(QP)Q 


( 6 ) 


where (Q ', P') is the image of (Q, P) under one inverse 
iteration of the map. U ( QP) is a function that depends 
only on the product QP and is analytic in a neighbor¬ 
hood of QP = 0. Thus the product QP is preserved 
by mapping of Eq. (6): Q'P' = QP which yields an in¬ 
variant hyperbola under the iterations, traced out by the 
dashed curves in Fig. 3. Under successive inverse itera¬ 
tions, A will follow the hyperbola outward to infinity, the 
distance between A' and the stable manifold will decrease 
by a factor of U(QP ) after each iteration, thus X will 
converge to the stable manifold exponentially. Transfer 
back to the phase space in (q,p) coordinates shows that 
this exponential convergence always holds [46]. 

Equivalently, it can be shown that any deviation inside 
the convergence zone of the unstable manifold will col¬ 
lapse onto the unstable manifold under forward iteration. 
The details are skipped here. This provides a scheme to 
characterize the behavior of points near the manifolds 
and implies the manifolds have a certain tolerance to ini¬ 
tial errors when generated numerically. This tolerance 
is the manifold stability which is the starting point of 
locating successive iterations of heteroclinic intersecting 
points. 


B. Heteroclinic orbits 

The method section is illustrated by calculating the 
orbit segment of the heteroclinic intersecting point 
i?o £ U+(0, 0) H 5 + (0.5,0) shown in Fig. 4. It facilitates 
a straightforward insertion technique to insert new 
points into the manifolds at each iteration to maintain 
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a dense enough set of points to better interpolate the 
manifolds. This kind of insertion technique (and more 
sophisticated ones) has been used by many authors to 
generate the invariant manifolds to a high accuracy 
on fine scales. See for example [29] for a widely used 
technique which sets up a criterion examining the local 
curvature of the manifold and inserts points if the 
curvature exceeds certain thresholds. A more elaborated 
method is developed in [47], which includes a point 
redistribution procedure after the insertion to achieve 
greater resolution in certain regions of interests. A 
notable generalization of these techniques into aperiodi- 
cally time-dependent vector fields is introduced in [39]. 
Approximation to the invariant manifolds using geodesic 
circles and the Hobson’s criterion in higher dimensional 
cases can be found in [36, 38]. Here, only the part of the 
manifolds near the homoclinic intersections are needed. 
Thus a linear insertion of new points at every iteration 
is sufficient, which is demonstrated in the following. 
Define the positive (negative) direction on C/ + (0,0) 
to be pointing outward (inward) from (towards) (0,0) 



FIG. 4. Heteroclinic tangle of the kicked rotor. Shown here 
are U + (0,0) and S + (0.5,0). Note that they are the same 
manifolds shown in Fig. 1, with C/ + (0,0) extending further. 
The phase point Ro in the lower zoomed-in graph is also the 
same point as in Fig. 1. 

along the manifold, and the positive (negative) direction 
of S + (0.5,0) to be pointing outward (inward) from 
(towards) (0.5,0) along the manifold. Calculate Rq 
first by interpolation methods. Then pick two points 
on U + (0, 0), namely and Uq 2 \ such that Uq 1 '* is 


( 2 ) 

adjacent to Rq on the negative side, and Uq is the point 
adjacent to Rq on the positive side. Another pair of 
points on S' + (0.5, 0) are picked in the same way, namely 

Sq 1 * and Sq 2 \ such that Sq ^ is the point adjacent to 

( 2 ) 

Ro on the negative side, and Sq is adjacent to Rq on 
the positive side. Due to the reflective property of point 
( 0 , 0 ), the Ri switch back and forth under each iteration 
as follows: 

R±i € C/ + (0,0) f) 5 + (0.5,0), (i even) 

R ±i e U~ (0,0) p| 5+(0.5,0), (i odd). 

The scheme is to obtain the set (i even) by in¬ 

tersecting Z7 + (0, 0) with A -1 " (0.5, 0) and iterate it once to 
obtain the set {i?±i} (i odd). 

Consider the inverse iterations with even iteration 
numbers first. R-i can be calculated by iterating Rq 
with the inverse map directly. However, the maximum 
number of iterations is limited due to the exponential 
growth rate of initial error associated with Rq. Under in¬ 
verse iterations, the component of the error transverse to 
the stable manifold will vanish by the virtue of manifold 
stability, but the tangential component will be magnified 
along the stable manifold. This will cause the calculated 
R-i to “drift away” from the true result exponentially 
fast. To avoid this, instead of iterating R 0: iterate Sq 
and Sq 2 ' 1 inversely to obtain their images S_] and S* 2 }. 
Then locate the intersecting point between f/ + (0, 0) and 
the straight line segment connecting S* 1 * and S (2 J. They 
provide a better estimate of R-i. 

There is a difficulty however, i.e. the exponential 
growth of the distance between S^J and S 1,2 *. An im¬ 
proved procedure is illustrated in Fig. 5. As the distance 
becomes larger, at some iteration the straight line con¬ 
necting the two points becomes a bad approximation to 
the local stable manifold (even though the propagated 
points, S^* and S* 2 *. continue to converge towards the 
stable manifold). Take the example of inverse iterations 
from Rq to i?_ 2 . Instead of iterating Sq X * and Sq"' 1 in¬ 
versely to obtain and S^.j , insert two new points 
5' (1) and Sq 2 ' 1 on the line segment connecting Sq ^ and 
Sq 2 \ such that the distance between Rq and Sq 1 * is re¬ 
duced by the instability 1 /U 2 {QP) relative to the dis¬ 
tance between R 0 and Sq 1 1 . Do similarly with Sq 2 ' 1 . Then 
iterate Sq ^ and Sq 2 * inversely twice to obtain a new pair 
of points S^ and ^- 2 - The structural stability ensures 

that any error placing Sq X * and Sq 2 * on the stable man¬ 
ifold collapses exponentially. 

The orbit point I ?_ 2 is located as the intersection be¬ 
tween the unstable manifold and the straight line con¬ 
necting S^ and S (2) 2 . Repeat the same insertion and 
iteration process to get R- 4 , R-e, Practical modifi¬ 
cations of the simple description of the process given here, 
such as using four or six iterations of the map instead of 
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FIG. 5. Calculating R-i -2 from R-i. Insert a new pair of 
points S'S^ and S'S^ near R-i, then iterate them inversely 

twice to <Sl_7_2 an( i S-i- 2 - R-i -2 is calculated as the inter¬ 
section between the unstable manifold and the line segment 
connecting 2 and S^_ 2 . The insertion of S'^ and S'^ 
is carried out in order to cancel out the exponential growth 
that otherwise would occur, thus ensuring a good approxima¬ 
tion to the local stable manifold. 

two or reducing the distance by more than 1/U 2 (QP), 
can be introduced for particular dynamical systems and 
circumstances. It is not necessary to invoke this entire 
procedure in order to obtain the odd iteration numbered 
points. One can just iterate R-i once inversely to ob¬ 
tain R_i_i. The accuracy of the heteroclinic orbit is not 
compromised in this way. 

The calculations of forward iterations are similar, the 
only difference is the use of unstable manifolds instead 
of stable ones. The idea is to find the pair of points Ly ' 1 
and U - 2 - 1 by the insertion of U'^\ and U'^\ on t/ + ( 0 , 0 ). 
The scheme is shown in Fig. 6 . This procedure calculates 
R+i ( i even), with the R+i ( i odd) calculated similarly 
as with the stable manifold by mapping the even iterates 
forward once. 


C. Kicked rotor heteroclinic orbits 

The above scheme applies to the heteroclinic orbit 
starting from Rq. By changing the initial choice of Rq 
within a single loop structure (or turnstile for the homo¬ 
clinic case) and repeated use of the scheme, we can cal¬ 
culate a whole set of distinct heteroclinic orbits. Orbits 
with longer excursion lengths only require that we gener¬ 
ate longer segments of the stable and unstable manifolds 
at the beginning, which do not introduce any difficulties. 
In this section numerical results are given for two example 
heteroclinic orbits. The first case is given in detail for the 
particular orbit Rq shown in Figs. 1,4, which cannot be 
followed directly with the mapping forward and inversely 
with time in a double precision calculation better than 
R -5 to R +,. Its initial condition is well approximated as 



FIG. 6. Calculating Ri +2 from R+. Insert a new pair of points 
U[^ and U near Ri, then iterate them forward twice to 
U}1 j and U -%. Ri +2 is calculated as the intersection be¬ 
tween the stable manifold and the line segment connecting 
[/}+ 2 an d Ri+ 2 - The insertion of and cancels out 

the exponential deviation between and U^ 2 \ ensuring a 
good approximation to the local unstable manifold. 

R 0 « (0.44217031018010239,0.51879785319959426). It 
happens to be the heteroclinic orbit with the simplest 
excursion away from {X a ,Xp\. With the method just 
described, it is possible without any additional techni¬ 
cal enhancements to follow this heteroclinic orbit from 
R_ i 9 to i?i 4 - The inverse and forward time limits are 
determined by the fact that the orbit expressed in dou¬ 
ble precision is indistinguishably close to X a and Xp by 
those times. As the orbit approaches fixed limiting points 
in its past and future exponentially quickly, it is not pos¬ 
sible to illustrate it very well in a simple phase space plot 
with successive points labeled {R -19 • • • R\a}. Instead, a 
schematic is shown of the first few iterations beginning 
with i?_ig in Fig. 7. 

To show that the solution is behaving properly, con¬ 
sider the stability analysis of the initial and final fixed 
points. In a small neighborhood of either X a = (0, 0) 
or Xp = (0.5,0), the mapping can be approximated 
by the linearized tangential equation 6 X a rp\(n + 1) = 
M XaW 6 X a (p)(n), where SX a (p)(n) = (Sq n , 6 p n )^ (/s) is 
the deviation of the n th iteration relative to X a rp\, and 
Mx a(l j) is the stability (Jacobi) matrix of X a (py A 
straightforward calculation of the stability matrices for 
both X a rp ) yields: 

6 X a (n + 1) = l) SXo(r ^ 

5Xp(n + l)= (^ + k K ^j 6 Xp(n) . (7) 

The eigenvalues associated with the stable manifolds of 
X a{ p) are \ s a{0) = (2 ± K =p \fK 2 ± 4 K) /2, where the 
lower sign refers to a and the upper sign to /3. Likewise, 
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In d n 
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-15 


-25 


-35 


8 12 16 
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FIG. 7. Distance between forward iterations of R-ig and 
(0,0) (not to scale). R~ig is marked by do and R- ig+„ by 
d n . Note that the integer n should be an odd number since 
the orbit point i?_ig +n is located on the upper branch in our 
figure. (7 + (0,0) is the solid line while R~(0,0) the dashed 
line. 


FIG. 8. Exponential convergence under inverse iteration near 
X a . The slope of the graph matches the characteristic expo¬ 
nent to five decimal places; ln|Aa| = 1.80593 versus 1.80592 
for the trajectory. 


the eigenvalues associated with the unstable manifolds 
are A“ (/3) = (2 ±K± sjK 2 ± AK)/2. Setting K = 8.25 
gives the four values relevant to this orbit example. 

For the early iterations R n , n close to —19, consider 
the norm of the difference between R n and X a , d n = 
\J5q2 + Sp^. It should turn out that 

d n +1 ~ |A“|d n (8) 

is an excellent approximation to the true dynamics of the 
heteroclinic orbit. Following the orbit for multiple itera¬ 
tions as in Fig. 7 and considering the natural logarithm 
of the norm leads to 

lnd„ ss rain |A“| + lnd 0 (9) 

Therefore a plot of In d n versus n should be nearly a 
straight line with slope ln|A“|. Shown in Fig. 8 is the 
graph plotted from our calculations of the orbit segment 
using the manifold intersections method. The points do 
indeed line up with the correct slope to a high degree of 
precision. 

Similar verification can be done for the forward itera¬ 
tions approaching Xp. Start with the point i?i 4 , which 
is the furthest forward iteration found, and count back¬ 
ward. Let do be the distance between and Xp, d n 
the distance between i?i 4 _ n and Xp. The same relation 
applies as before except using |A^|. The plot looks iden¬ 
tical to that of Fig. 8, so is not included here. The only 
difference is that it has a slope of 2.31808, which is to be 
compared with In |Ao| = 2.31762. This also differs in the 
fifth significant digit. 

As mentioned in the introduction, in semiclassical the¬ 
ory classical action differences divided by H control quan¬ 
tum interferences. Thus, even very small errors in the 
calculation of action differences, depending on the value 


of H, can ruin the quality of a semiclassical approxima¬ 
tion. Exploiting the exact relationship relating the ac¬ 
tion difference between two heteroclinic orbits to the ge¬ 
ometric area under the invariant manifolds that connect 
them [48, 49] gives a meaningful measure of the quality 
of an orbit’s calculation. 

For the kicked rotor, the action function can be written 
explicitly as: 


F(X n , X n+1 ) = {<ln+1 2 gn)2 + cos 27 rq n (10) 

where X n = (q n ,Pn) is an arbitrary phase point that 
is mapped to some X n+ i = (q n+ i,p n+ i) under Eq. (4) 
without the mod operation. If used as an Fi(g, Q\ t) gen¬ 
erating function, — F generates such a map. The two 
classical action functions F{X a ^p\,X a ip^) = ±I\/Air 2 are 
invariant under translation in time as they refer to fixed 
points of the map. Repeated use of Eq. (5.6) from [49] 
gives: 

= J Pdq ( 11 ) 

U[X a ,R 0 ] 


E 


F(Ri-i,Ri) ~ F(X a , X a ) 


The path of integration U[X a ,Ro] is the segment of un¬ 
stable manifold from X a to Rq, as shown in Fig. 9. The 
direction is denoted by an arrow. Similarly: 


+oo 

E 


2—0 


F(Ri,R i+ 1 )-F(Xp,Xp) 



S[R 0 ,X P ] 


( 12 ) 


The path of integration S[Ro , Xp] is the segment of stable 
manifold from Rq to Xp, as shown in Fig. 9. Combining 
























0.6 


1 



q 


FIG. 9. Path of integration to calculate the geometric ar¬ 
eas. U[X a ,Ro] is the segment of the unstable manifold from 
X a to J?(0). 5[I?o, Xp] is the segment of the stable manifold 
from Ro to Xp. The shaded region marked ‘A’ shows the 
corresponding area. 


the above two equations gives: 


E 


F(Ri-\,Ri) - F(X a ,X a ) 


+ 


1 =— oo 


+oo 

E 


2 = 0 


F(Ri,Ri+i)-F( Xp,Xp) 



U[X a ,R 0 \ 


(13) 


For the classical action calculation, the infinite sums 
are truncated to R-ig and Ru, which are as close 
to the fixed points as a straightforward double preci¬ 
sion calculation can be. The classical action difference 
gives 0.12938887802084850, whereas a construction of 
the unstable and stable manifold segments and numer¬ 
ical integration of the area denoted ‘A’ in Fig. 9 gives 
0.12938887802085794. They differ in the \d th decimal 
place, the limit of double precision calculations. 

A more complicated heteroclinic orbit is denoted with 
the point Rq shown in Fig. 10. It not only makes a 
more complicated excursion away from {X ai Xp}, it also 
wraps around both directions of the unfolded torus, and 
thus translates from one phase space copy to another. 
The peculiarity of this kind of heteroclinic orbit is that 
its fixed point Xp with nonzero rip-winding number is 
constantly shifting in the q coordinate on the phase 
plane tiled with unfolded tori under iteration, along 
with its stable and unstable manifolds when the map 
is applied. This follows from the mapping equations, 
Eq. (4), by dropping the mod operation. In this case, 
the point (—1.5, —2) is mapped to (—1.5 — 2n, —2) under 
n iterations. The stable and unstable manifolds of it are 
shifted the same way. Thus: 


0 



-2 


-3 




FIG. 10. Heteroclinic intersection Ro between 

F + (0,0) and S + (— 1.5, — 2). Note that S + (— 1.5, — 2) 
is just S + (0.5,0) (same Xp) shifted by the wind¬ 
ing numbers (-2,-2). Its initial values are Ro = 
(-1.6105430949740283,-1.0599500106337416). 


-Ro £ S' + (—1.5, —2) P| U(0, 0), 

R n e S+(- 1.5 - 2n, -2) f| 1/(0,0) 

It was possible to obtain inverse iterations up to R-n 
and forward iterations up to R 13 . Similar to the previous 
case, the values of In d n versus n compared quite closely 
with the characteristic exponents of X a and Xp. The 
beginning slope from the values of In d n is 1.80589, which 
is the same to 5 decimal places with ln|A“| = 1.80594. 
Likewise, but a bit less accurately, the later slope from 
the values of In d n is 2.32168, which is roughly the same 
to 4 decimal places with In |A^| = 2.31762. 

In order to study the area-action relation for this or¬ 
bit, a slight modification of the algorithm is needed to 
account properly for the shifting of the phase space un¬ 
folded torus. Let the path U[X a ,R 0 ] be the segment of 
l/+(0, 0) from X a to i?(0), and the path S[Rq, Xp] be the 
segment of S + {— 1.5, — 2) from R 0 to Xp = (—1.5,—2). 
For the U[X a , i?o] path of the heteroclinic orbit, there is 
no definitional change from the previous case, just the use 
of a different orbit; i.e. F(i?,_ 1 , Rt) is the action function 
evaluated from Ri-\ to Ri and F(X a , X a ) is the action 
function of the fixed point X a , independent of i and equal 
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to K/Ak 2 . 

For the latter half of the path, i.e. let 

F(Xp, Xp) be the action function that maps Xp = 
(-1.5 - 2i, —2) to Xp = (-1.5 - 2(i + 1), —2) as the 
map takes the point Ri to R.-j+x- The modified action is 

F(Xp,Xp)=2-^,\/i. (14) 

The sum of the actions that contribute to double preci¬ 
sion are: 


E 

2 = —16 


F{Ri-\,Ri) — F(X a ,X a ) 

+ 12 


+ 


E 

i=0 


F(R tl R l+1 )-F(Xp,Xp) 

= 0.16465128640816951 



q 


FIG. 11. Relevant areas Ai, A 2 and A 3 for the action cal¬ 
culation. U[X a ,Ro] is the segment of t/ + (0,0) from X a to 
R 0 - is the segment of S + (— 1.5, — 2) from R 0 to 

Xp. The path integration f pdq is the alge- 

U[X a ,fl 0 ]+S[«o,^,3] 

braic area —A\ + A 2 — A 3 . 

The geometric area for this orbit is less straightfor¬ 
ward to conceptualize. It is shown in Fig. 11 for clar¬ 
ity. Construction of the stable and unstable manifolds 
and using a numerical integration of the relevant areas of 
—A 1 + A 2 — A 3 gives a total of: 

J pdq = - A 1 + A 2 - A 3 

U[X a ,Ro]+S\R 0 ,Xfi] 

= 0.16465128641878213 

Deviation of the two numbers begins in the 11 th decimal 
place. This is not quite as accurate as the first hete¬ 
roclinic orbit case, the simpler one. However, the area 
integral is more difficult to calculate accurately and it 
is not known which source of possible error, either from 
the orbit action sum or the numerical area calculation, 
gives rise to the increased inaccuracy. Nevertheless, the 
accuracy is excellent. 


D. Action relations between fixed points 


Here a general formula relating the action difference 
between a pair of fixed points to a region bounded by 
certain segments of the unstable and stable manifolds is 
derived. The result obtained here is independent of the 
context of the kicked rotor and applies to any Hamilto¬ 
nian system with a heteroclinic tangle. 

Equation (13) corresponds to the case that the switch¬ 
ing point from the unstable to the stable manifold along 
the integration path is chosen to be Rq. In general, it 
is possible to change the switching point to other orbit 
points Rk from {i?o} ; and modify the integration paths 
correspondingly. Thus a more general relationship can 
be given: 



F(Ri~i, Ri) — F(X a , X a ) 


+ 


+ OO 

E 


i—k 


F{R u R l+l )-F{Xp,Xp) 



U[X a ,R k \ S[R k ,X fi ] 


(15) 


Calculating the difference between Eq. (13) and 
Eq. (15) gives the action difference between the two fixed 
points X a and Xp. Subtracting Eq. (13) from Eq. (15) 
leads to: 


k ■ 


F(Xp, Xp)—F(X a , X a ) 


J Tdg - J pdq 

U[X a ,R k ] U[X a ,R 0 ] 

+ / / 

S[R k ,Xp] S[R 0 ,Xp] 


(16) 


Pdq 


For the cases that both X a and Xp are non-reflective, 
U[X a ,Rk] and U[X a ,R 0 ] are segments belonging to the 
same branch of the unstable manifold of X a , therefore: 



U[X a ,R k ) C/[X„,fio] 

and similarly: 



S[R k ,Xp\ S[Ro,Xp] S[R k ,Ro] 

Equation (16) simplifies to: 


k- 


F(Xp,Xp)-F(X a ,X a ) 



U[R 0 ,R k ] S[R k ,R 0 ] 


( 19 ) 
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which holds true for Hamiltonian systems with non- 
reflective fixed points. 

In the presence of reflective fixed points, such as the 
kicked rotor, U[X a ,R k ] and U[X ai R 0 ] are on the same 
branch only when k is an even number (twice iterated 
map). Letting k = 2, we have: 


J pdq+ J pdq = —B 
U[Ro,R2] S[R 2 ,Ro) 


= 2 


F(X p ,X p )-F(X a ,X a ) 


( 20 ) 


K 

-77-2 


where the B is the area of the phase space region en¬ 
closed by U[R 0 ,R 2 ] and *5*[7? 2 , -Ro] shown in Fig. 12. No¬ 
tice that the union of [/[-R 07 -R 2 ] and S'[f? 2 ,^o] gives the 
fundamental loop structure of the heteroclinic tangle 
under the twice iterated map. This is in the sense that 
the loop is the smallest object which can be used to 



q 

FIG. 12. The phase space region enclosed by U[Rq,R 2 ] and 
S[i? 2 ,-Ro] is labeled by B. Since U[Rq,R 2 \ and 5'[l?2,l?o] are 
the fundamental segments of the twice iterated map, in this 
heteroclinic tangle, B plays the role of the turnstile structure 
commonly seen in homoclinic tangles. Note that R 2 is ex¬ 
tremely close to Xp, making it hard to distinguish between 
the two in the scale of the figure. 


generally. Numerical calculation of the area gives: 

— B = -0.83589976498504137 (21) 

and 

- — = -0.83589976504928665 (22) 

7 T z 

They match up to the 11 th decimal place, which verifies 
Eq. (20). 


E. Heteroclinic fundamental loop structure 

There is necessarily a marked difference between the 
heteroclinic fundamental loop structure and the homo¬ 
clinic fundamental loop structure (the turnstile). It is 



be mapped forward and inversely in order to generate 
the full heteroclinic tangle, L/ + (0,0) and S ,+ (0.5, 0), and 
on the loop each distinct heteroclinic orbit has only one 
intersection point. In the common case of homoclinic 
tangles, the fundamental loop form a turnstile structure 
which was extensively used to study the flux entering and 
leaving certain phase space regions [48, 50, 51]. However, 
in our case the turnstile structure is replaced by the loop 
shown in Fig. 12, which is topologically equivalent to a 
circle. The area of the loop is equal to twice the ac¬ 
tion difference between the fixed points. Thus, the loop 
structure of the heteroclinic tangle gives the action dif¬ 
ference (or multiple thereof) of two periodic orbits more 


FIG. 13. The top panel shows the turnstile of the homo- 
clinic tangle formed by U + (— 0.5, 0), S + (0.5,0) ,S~ (—0.5,0) 
and U~ (0.5,0), whose two lobe areas must cancel. Structures 
like this are closely linked to the transport between differ¬ 
ent regions of phase space and are frequently mentioned in 
research literature. The bottom panel, which shows the fun¬ 
damental loop of the heteroclinic tangle between (0, 0) and 
(0.5,0), is quite different. It doesn’t form a turnstile shape 
and is equivalent to a single circle. 

well-known that the fundamental loop of the homoclinic 
tangle must have at least one crossing point between the 
stable and unstable manifolds in order to cancel out the 
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flux in and out of a certain phase region [48, 50, 51]. Us¬ 
ing the action relations between the fixed points of Sub¬ 
section HID implies the same conclusion. Furthermore, 
for the case of the heteroclinic tangle, the fundamental 
loop doesn’t even need a crossing point between the sta¬ 
ble and unstable manifolds, as shown below. 

For both the case without reflection (Eq. (19)) and the 
case with reflection (Eq. (20)), the action difference be¬ 
tween two fixed points is expressed as the algebraic area 
of the fundamental loop structures. Therefore for ho¬ 
moclinic tangles in which the two fixed points coincide, 
the algebraic area of the fundamental loop is necessarily 
zero, ensuring zero difference in the actions of the fixed 
points. This leads to the “8” shaped turnstile structure 
in the simplest case, or structures with more crossings 
between the stable and unstable manifolds in order to 
maintain a zero algebraic area. For heteroclinic tangles 
in which the two fixed points have different action func¬ 
tions, the algebraic area of the fundamental loop cannot 
vanish. For the kicked rotor, this circle shaped structure 
is already demonstrated by Fig. 12. 

The topology of the loop in Fig. 12 can be better il¬ 
lustrated if we shorten the unstable segment by iterating 
the loop inversely for one iteration, as shown in the lower 
panel of Fig. 13. The length of the stable and unstable 
segments are now closer to each other after the iteration, 
and the single loop is equivalent to a circle. The upper 
panel of Fig. 13 is the well known turnstile structure of 
the homoclinic tangle formed by the upper branches of 
the unstable manifold of (—0.5,0) and stable manifold of 
(0.5, 0), which is in clear contrast to the lower panel. 

IV. CONCLUSIONS 

There are many contexts in which one might be inter¬ 
ested in the actions of periodic, heteroclinic, or homo¬ 
clinic orbits. It is worth noting that in the context of the 
trace formulae or semiclassical propagation of wave pack¬ 
ets for chaotic systems, accurate calculation of periodic, 
heteroclinic or similarly homoclinic orbits is necessary for 
getting delicate quantum interference phenomena right. 
A simple method to calculate such orbits in a strongly 
chaotic system is given here. It works far better than 
any attempt to follow such a trajectory forward and in¬ 
versely in time directly using the equations of motion. 
The method relies on the structural stability of the in¬ 
variant stable and unstable manifolds, and the exponen¬ 
tial convergence of phase points in their neighborhoods 
towards them (unstable manifold forward in time, stable 


manifold inversely in time). The high degree of accuracy 
of the method was demonstrated with the help of a stabil¬ 
ity analysis of the limiting phase points ( X a , Xp) as well 
as an exact relation between certain classical action dif¬ 
ferences and their equivalence to the related phase space 
areas. 

An extremely interesting relation follows from the con¬ 
nection between phase space areas and classical action 
differences. For the heteroclinic case, assuming a generic 
case in which the two fixed points have different actions 
as in the example shown, the closed path along the fun¬ 
damental heteroclinic loop structure must enclose a non¬ 
vanishing area. Remarkably, the intersections of the sta¬ 
ble and unstable manifold has periodic orbit action differ¬ 
ences built into it. Contrast this to the homoclinic case 
where the loop structure is a turnstile. There the closed 
path winds around in a figure eight path and the total 
enclosed area must vanish because X a and Xp belong to 
the same orbit and their action difference vanishes. A 
heteroclinic loop structure cannot be a turnstile whose 
closed loop path integral vanishes. 

The demonstration was carried out for a paradigm of 
chaos studies, the kicked rotor. It takes the form of a dis¬ 
crete time dynamical map, but the same method could 
equally well be applied to continuous time dynamical 
systems by just creating an appropriate Poincare map. 
Higher dimensional generalizations of the method is di¬ 
rect by using more than two points each for the inter¬ 
polation of local unstable and stable manifolds. For ex¬ 
ample, in the case of a two-dimensional stable manifold 
intersecting a two-dimensional unstable manifold in four- 
dimensional phase space, the local shape of the manifolds 
are interpolated by four points on the tangent plane, so 
intersection can be located by intersecting the tangent 
planes. Better accuracy is obtained by using denser set 
of points on the manifolds or more sophisticated polyno¬ 
mial interpolation techniques. The applicability is there¬ 
fore rather general. 
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